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We derive the non-Maxwellian distribution of self-gravitating iV-body systems around 
the core by a model based on the random process with the additive and the multiplica- 
tive noise. The number density can be obtained through the steady state solution of 
the Fokker-Planck equation corresponding to the random process. We exhibit that the 
number density becomes equal to that of the King model around the core by adjusting 
the friction coefficient and the intensity of the multiplicative noise. We also show that 
our model can be applied in the system which has a heavier particle. Moreover, we 
confirm the validity of our model by comparing with our numerical simulation. 
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How relevant is the equilibrium statistical mechanics when we describe the steady 
state of a self-gravitating system (SGS) where many particles interact via the grav- 
itational force? Let's assume that the state of the SGS with equal mass m becomes 
isothermal with temperature T and the particles of the system will be distributed 
spherically symmetrically. Then, the structure in the phase space can be determined 
by the Maxwell-Boltzmann distribution. For example, the number density at a radial 
distance r in the real space is 

n MB (r) ex e~^* (r) , (1) 

where $(r) is the mean gravitational potential per mass generated by this whole sys- 
tem at r and k-& is the Boltzmann constant. This potential should satisfy a relation 
with the number density by the Poisson equation A$(r) = 4iTGmnMB{f) where G is 
the gravitational constant. A special solution of eq.(l) and this Poisson equation is 
^mb(^) = ksT /2nGm 2 r 2 known as the singular isothermal sphere. ^ This solution has 
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two problems: infinite density at r = and infinite total mass. Even though we solve 
the equations with finite density at r = 0, the solutions behave oc r~ 2 at large r, and so 
we cannot get around the infinite total mass problem. In either case, the solutions are 
unrealistic. As expected, it is unreasonable to apply the standard statistical mechanics 
to SGS where the interaction is a long-range force. 

However, the real systems in the universe, e.g. globular clusters, galaxies etc., have 
various structures. As for the most of globular clusters, it is known that the number 
densities of them in the real space have a flat core and behave as a power law outside the 
core. King interpreted these profiles by introducing the new distribution function (DF) 
in the phase space, known as the lowered Maxwellian, which becomes zero when the 
total energy is greater than a certain value by subtracting a constant from the original 
Maxwell-Boltzmann distribution. This is called as the King model. 2 ) The approximation 
of the number density of the King model, n,KM( r ), around the core is 

n KM(r) oc (1+r2 1 /fl2)3/2 , (2) 

where a is the core radius. Since he put forward this model, the number density has 
been applied to fitting for the surface brightness of many globular clusters, for example 
as in ref. 3 -* 

In this Letter, we try to derive this non-Maxwellian distribution in the real space 
around the core from a new point of view. Our simple model uses a Brownian dynamics 
described by random process with the additive and the multiplicative noise, which is 
quite different from the King model because his procedure was done to the DF in the 
steady state. Physically, this model represents that the gravitational force induced by 
<3>(r), the mean force, fluctuates in time. We get the Fokker-Planck equation from the 
Langevin equation and show that the same result as the King model can be obtained 
from the steady state solution. In addition, we show that our model can deal with the 
case that the system includes another particle whose mass is M(> to), corresponding 
to the black hole in a globular cluster. 

Now, we investigate the steady number density (SND) of the SGS with mass to 
including a particle with mass M by using iV-body simulation. Especially, we show the 
results that the density profiles in the steady state have a core and behave as a power 
law. The system is composed of iV = 10000 particles. At t — 0, all velocities of the 
particles are zero and they are distributed by n (r) oc (1 + r 2 /a p 2 )~ 5 / 2 (0 < r < Aa p ) 
which is the density in the real space of Plummer's solution. ^ In this SGS we put 
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another particle with mass M in the origin at t — 0. We shall change the mass as 
M/m = 1, 5, and 10. Throughout this Letter, we adopt a unit system where the core 
radius of the Plummer's solution a p , initial free fall time tff, and the total mass N • m 
are unity. 

We started the iV-body simulation under these conditions. For dynamical evolution, 
we use GRAPE-7, special purpose computer for gravitational force. 4 * 1 For computation 
of the gravitational force, we apply Plummer's softening: the potential energy between 
the ith and the jth particle separated by a distance is —Gm 2 / ^r^ 2 + e s 2 where e s is 
the softening parameter. We set e s = 1CT 3 . For evolution, we use sixth order symplectic 
integrator. 5 ^ The time step for the simulations is defined as At = 1CT 5 . We carried out 
the simulations until t = 100 tff. During simulations, the error of the total energy is 
less than 0.1%. 

At first, most of the particles collapse into the origin within several tff. Approxi- 
mately after 20 tff, the distribution becomes stable and the system goes to the steady 
state. In Fig.l, we show the logarithm of SND as a function of logr for M/m = 1, 5, 
and 10. For each M, the SND has a core and behaves as a power low at larger r than 
the core radius. 

We now fit SNDs around the core by ra fit (r) = C/(l + r 2 /a 2 ) /3 . The results are 
summarized in Table I. For M/m = 1 and 5, (5 ~ 3/2 which is similar to the exponent 
of the King model. The density at the origin C increases as M is increased, which is 
simply understood as a result that many particles are attracted by the heavier particle. 

Table I. The best fitting parameters of a function (7/(1 + r 2 /a 2 ) /3 for steady number densities shown 
in Fig.l 



M/m 


a 
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6.44 x 10~ 2 


1.49 


6.91 x 10 5 
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6.28 x 10~ 2 


1.49 


7.41 x 10 5 


10 


5.84 x 10~ 2 


1.45 


7.99 x 10 5 



In order to explain these results and derive the non-Maxwellian distribution around 
the core, we demonstrate the simple model based on Brownian motion which is quite 
different with the King model. The reason why Brownian motion appears in the SGS 
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Fig. 1. (Color online) Logarithm of the steady number density n(r) as a function of logr for (a) 
M/m — 1, (b) M/m = 5, and (c) M/m — 10. In each figure, the dashed curve and symbols denote a 
fitting curve n^(r) = C/(l + r 2 /a 2 ) 13 and the result of our numerical simulation, respectively. 



is as follows. After the collapse, the density around the core becomes high. Thus the 
particles around there disturb the orbits of others repeatedly, so that the movements 
become random. 6 ) As the time to occur this disturbance, we introduce the local two- 
body relaxation time £ r ei' 7 ^ 

0.065(i(r) 3 

t Ie \(r) = — , (3) 

G 2 n(r)m 2 ln(l/e s ) 

where a(r) is the standard deviation of the velocity at r and we adopted ln(l/e s ) as 
the Coulomb logarithm. 

Figure 2 shows the logarithm of t re i, which is calculated by using er(r) and n(r) 
obtained from our iV-body simulation, as a function of logr. As expected, t re \ around 
the core is short. Our simulation continues after the collapse during about 80 tff which 
is sufficiently longer than t Te \ around the core. As the radius increases, however, t ie i 
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becomes longer than the rest of our simulation time, which means that the Brownian 
motion does not occur at large r. Therefore, note that our model is valid only in the 
neighborhood of the core. 
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Fig. 2. (Color online) Logarithm of the local two-body relaxation time t re \ as a function of logr for 
M/m = 1, 5, and 10. 

When constructing our model, the following points are premised: the model describes 
the dynamics near the steady state and the mean distribution is spherically symmetric. 
As well known, a gravitational force at r arising from such a spherically symmetric 
system depends only on the particles which exist inside the sphere with radius r and this 
attractive force acts along the radial direction. In other words, this mean force — /(r) is 
the gradient of the mean potential: —f(r) = —md<&(r)/dr. Of course, lim r _> /( r ) = 0. 
Hence, we assume 

f{r) = am 2 r + 0{r 2 ) , (4) 

where a = 47rC7C/3. This condition is necessary in order that the SND has a core as 
will become clear later. Conversely, if the SND has a core, f(r) is represented as eq.(4), 
which is easily understood as following way. Now, the SND is almost constant in the 
neighborhood of the core, and so the sum of the mass of particles which exist inside 
small r is M(r) ~ 4nCmr 3 /3. Therefore, f(r) = GmM(r)/r 2 ~ 4irGCm 2 r /3. After all, 
the SND is linked to the mean force f(r) self-consistently around the core. 

For the case M/m = 1, we can identify another particle with others. Contrary to 
this, we must consider the effect of the particle for the case M/m ^ 1. Now, we suppose 
that the heavier particle exists at the origin. Then, the attractive force by this particle 
at r is —F(r) = —GmM/r 2 . We can estimate f(r) around there as f(r) ~ 4TrGm 2 Cr/3. 
Thus, F(r)/f(r) ~ 3Mr~ 3 /AwCm. This ratio has a meaning when r ~ 10~ 2 . Therefore, 
if r is smaller than the radius, particles are influenced by not only f(r) but also F(r), so 
that the core should disappear. In fact, we have done a numerical simulation with the 
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heavier particle fixed at the origin, where this result is confirmed. On the other hand, 
Miocchi improved the King model in order to describe steady state of a globular cluster 
including an intermediate-mass black hole and reported that the density becomes cuspy 
as the mass of the black hole increases. 8 -* Because the heavier particle of our numerical 
simulation is not heavy so much, the particle is not trapped at the origin. Therefore, 
we do not consider the effect of the heavier particle explicitly and we suppose that the 
particle influence SGS through the density at the origin C: as M becomes larger, it 
attracts more particles and C increases as shown in Table I. Thus, a = AnGC/3 is 
increasing function of M. 

It is natural to consider that the distribution fluctuates around the mean value 
because of the many disturbances. The fluctuating part of distribution should not be 
spherically symmetric, so that this produces the forces along not only the radial di- 
rection, but also the other directions. We assume that they are random forces and set 
the intensity of them at r 2E(r) 2 . In addition to such random forces resulting from the 
fluctuating distribution, a particle at r is expected to be influenced by random forces 
generated from neighbor particles. We set the intensity 2D which is independent with 
position. 

Brownian motion with the above assumptions is described by the following Langevin 
equations in spherical coordinates: the radial direction 

m 7 r = -f(r) + V2E(r)r] r (t) + V2D£ r (t) , (5) 

the elevation direction 

mi r6 = V2E(r)r]e(t) + V2D&(t) , (6) 

and the azimuth direction 

m 7 r sin 6<j) = v / 2 J B(r)^(t) + V2D^(t) , (7) 

where 7 is the coefficient of dynamical friction in the low velocity limit and independent 
on the velocity. 1 - ) In the Chandrasekhar dynamical friction formula, the coefficient is 
more complicated. 1,9 ) But we use the coefficient in such a limit, because the density 
around the core is so high that particles around there move slowly. 

Since we consider the Brownian dynamics near the steady state, the inertial terms 
are neglected. 10 ) The noises in each Langevin equation, £ and r/, are zero- mean white 
Gaussian and are correlated only to themselves. Of course, the correlation function is 
the Dirac delta function. 12 ) 
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Here, we make the following final assumption: the intensity E(r) behaves similarly 
to /(r) around the core, that is E(r) = y/ef(r) = y^eam 2 r + 0(r 2 ) where e is a positive 
constant. Then, the first and the second term on the right-hand side of eq.(5) become 

-f(r){l - V2~e Vr (t)} . (8) 

This represents that the mean force fluctuates. As we wrote before, the distribution in 
the near steady state is expected to fluctuate around the mean value which yields the 
mean force —f(r), so that the total gravitational force from the whole system along the 
radial direction is described as in eq.(8). 

We have the Fokker-Planck equation governing the spherically symmetric probabil- 
ity distribution function (PDF) P(r, t) 

—P(r t) - {— + - — \ P(r t)-\ — — — r 2 f(r)P(r t) 

dt ' (m / y) 2 \ dr 2 r dr J ' m^y r 2 dr ' 

(9) 

where the prime indicates derivative with respect to r. Then, the PDF with the Jacobian 

p(r,t) = 4nr 2 P(r,t) satisfies the following Fokker-Planck equation. 

d , v D f d 2 8 2}., 1 w x 

dt P ^ t] = (^l^-^r/ /?(r ' t) + ^^ /(r)/?(r ' t) 

+ (^{l^^-|^ 2 }^^ < 10 > 
This equation is useful when integrating with respect to r. 

The steady state solution p s t( r ) satisfies the equation (10) with the left-hand side 
zero. By integrating the equation with respect to r, we have 

{(^ + (^ /(r) 1 /4(r) 

' f(r)f'(r)--f(r) 2( 



Pst(r) 



_{m^j) 2 r (rrry) 2 { r J rwy 

= const. . (11) 

Now, we impose the binary condition that P st (r) = p st (r)/(47rr 2 ) and the derivative 
do not diverge at the origin. Then, when r — > 0, p st (r) = 0(r 2 ) and lim r ^ Pst( r ) = 
lim^o 47r(2r P st (r) + r 2 P s ' t (r)) = by which the constant on the right-hand side of 
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eq.(ll) is decided and we obtain 



d (r) - rf\r){ef\r)+ mi }-2{D + ef(rf} 

P«lr) ~ r{D + efm Pst(r) • (12) 



Let's represent f(r) by the first order of r. Then, 



Here, if we set 



' ( r \ — ea 2 m 4 I ema / 

stK')- f D 

2_ ^ _ J Q _ 1 ^ 7 



Pst(r) . (13) 



a 



and /3 = 1 (JL- + l) , (14) 



the equation (13) can be solved like 



r 2 



« ' (15) 



which yields 



PstW K (1 + rVaT ' (16) 
Since the stochastic behavior of a particle of the SGS is governed by this PDF, the SND 

of this system can be obtained by multiplying this equation by a constant. 

As in eq.(14), the exponent (3 must be larger than 1/2, which does not contradict 
our numerical simulation shown in Table I. In order that the equation (16) corresponds 
completely to the King model, /3 = 3/2or7 = 2ema must hold. We can regard this 
relation between the friction coefficient 7 and the intensity of the multiplicative noise 
e as a kind of fluctuation- dissipation relation, 13 ^ which usually plays an important role 
when a random process with an additive noise goes to the equilibrium state described 
by the Maxwell-Boltzmann distribution. 

The core radius a is proportional to a square root of the intensity of the additive 
noise D owing to eq.(14). Then, the intensity spreads the region where the density is 
almost constant. This is recognized as the effect of this noise which has the property 
that it makes a system homogeneous. 

Now, let's examine the role of the mass of the heavier particle M in this system by 
the naive discussion. As written previously, a is increasing function of M. a exists in the 
denominator of a and j3. Then, both values should be reduced when M is increased if 
other parameters are independent on M. These theoretical expectations are consistent 
with our numerical results shown in Table I. 

How the result eq.(16) changes if the mean force does not fluctuate? The steady 
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state solution of eq.(9) with e = is 

P st (r) oc e -^ $ W . (17) 

Therefore, our result goes to the singular isothermal sphere as discussed at the beginning 
of this Letter when the mean force does not fluctuate. 

Here, we examine the relation between the King model and our model. King trans- 
formed the DF in the phase space in order to avoid the singular isothermal sphere. 
In our model, we introduce the multiplicative noise into the system influenced by the 
mean force and the additive noise whose PDF becomes Maxwellian in the steady state 
as shown in eq.(17), so that the non-Maxwellian DF eq.(16) is derived. In short, al- 
though these procedures seem to be different, they may have the same meaning at least 
around the core. 

In conclusion, we have derived the non-Maxwellian DF eq.(16) by the Brownian 
dynamics with the fluctuating mean force and the additive white noise. The number 
density derived from the DF can be the same as that of the King model around the 
core by controlling the friction coefficient and the intensity of the multiplicative noise. 
Furthermore, our model can be valid in the SGS with a heavier particle. Of course, these 
results are consistent with our numerical simulation. We can say that such a stochastic 
dynamics occurs behind the background of the King model. Finally, note that our result 
is available only in the neighborhood of the core. So, we must derive the density globally 
by further extended model and investigate the difference between the model and the 
King model, in the near future. 
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